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Abstract 

In this paper, we present a framework based on the generahzed lattice-Boltzmann equation 
(GLBE) using multiple relaxation times with forcing term for eddy capturing simulation of wall 
bounded turbulent flows. Due to its flexibility in using disparate relaxation times, the GLBE is 
well suited to maintaining numerical stability on coarser grids and in obtaining improved solu- 
tion fidelity of near-wall turbulent fluctuations. The subgrid scale (SGS) turbulence effects are 
represented by the standard Smagorinsky eddy-viscosity model, which is modified by using the 
van Driest wall-damping function to account for reduction of turbulent length scales near walls. 
In order to be able to simulate a wider class of problems, we introduce forcing terms, which can 
represent the effects of general non- uniform forms of forces, in the natural moment space of the 
GLBE. Expressions for the strain rate tensor used in the SGS model are derived in terms of the 
non-equilibrium moments of the GLBE to include such forcing terms, which comprise a general- 
ization of those presented in a recent work (Yu et al, Comput. Fluids, 35, 957 (2006)). Variable 
resolutions are introduced into this extended GLBE framework through a conservative multiblock 
approach. The approach, whose optimized implementation is also discussed, is assessed for two 
canonical flow problems bounded by walls, viz., fully-developed turbulent channel flow at a shear 
or friction Reynolds number (Re) of 183.6 based on the channel half-width and three-dimensional 
(3D) shear-driven flows in a cubical cavity at a Re of 12,000 based on the side length of the cavity. 
Comparisons of detailed computed near- wall turbulent flow structure, given in terms of various 
turbulence statistics, with available data, including those from direct numerical simulations (DNS) 
and experiments showed good agreement. The GLBE approach also exhibited markedly better sta- 
bility characteristics and avoided spurious near-wall turbulent fluctuations on coarser grids when 
compared with the single-relaxation-time (SRT)-based approach. Moreover, its implementation 
showed excellent parallel scalability on a large parallel cluster with over a thousand processors. 
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I. INTRODUCTION 



The lattice Boltzmann method (LBM), employing minimal discrete kinetic models to 
solve fluid mechanics and other physical problems, has attracted much attention in recent 



years 



mm. 



Instead of directly solving the Navier-Stokes equations (NSE), the LBM 
involves the solution of the lattice Boltzmann equation (LBE) 5|, la, [Zl, ISl, l9| , which describes 
the evolution of the distribution of particle populations on a lattice whose collective be- 
havior asymptotically reproduces the dynamics of fluid flow. More speciflcally, the lattice, 
possessing sufficient rotational and other symmetries, restricts the collisions and movements 
of particle populations along discrete directions, as represented by the LBE, in such a way 
that in the continuum limit, fluid flow represented by weakly compressible NSE is recovered. 
While its origins lie in the lattice gas cellular automata (LGCA) [10], its formal connection 



to kinetic theory ll|, [12] has more recently led to improved physical modeling using the 



13j, and greater amenability for numerical 



LBM, for example to represent multiphase flows 
analysis The attractiveness of the LBM comes from the simplicity of the stream-and- 
collide computational procedure, absence of the need for an elliptic Poisson-type equation 
for the pressure fleld, ease in handling boundary conditions for representation of complex 
geometries, and excellence parallel performance due to its explicit and local nature. As a 

nnn 

result, it has found a number of interesting fluid flow applications |2|. I4J. |15I|. 

Representation and computation of turbulence is one of the most challenging aspects 



of fluid dynamics la, ll7|]. In recent years, signiflcant progress has been made to derive 



turbulence models a priori from discrete kinetic theory 



18|, [l9|, 120] , and turbulence modeling 



in the LBM has found much success in practical applications, for e.g., by Teixeira 



2l| and 



Chen et al. [22]. Also, various prior studies have found that LBM is a reliable and accurate 



method for direct numerical 
- see for e.g. Refs. 



23 



24 



simulation (DNS) of various benchmark turbulent flow problems 



25 



27 



28 



30|. 



On the one hand, turbulence models in the Reynolds-averaged contexts are generally 
required to represent physics over a wide range of scales. While turbulence at small scales 
tends to be somewhat more universal, large scale turbulent motions are strongly problem 
dependent. Hence, it is unrealistic to expect Reynolds- averaged models to accommodate 
and represent large-scale behavior of different classes of turbulent flows in the same manner 
without resorting to considerable empiricism. On the other hand, the DNS approach resolves 
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all relevant spatial and temporal scales and can thus predict all possible fluid motions with 
high fidelity. However, its computational cost limits its utility to low Reynolds numbers. 
Thus, it is often more practical to use large eddy simulations (LES), where fluid motions with 
length scales greater than the grid size are computed and the effect of the unresolved eddies 
at subgrid scales (SGS) are modeled [Sli- In this regard, while the use of simple Smagorinsky 
model {32! to represent SGS effects and perform LES using LBM was proposed some time 



ago by Hou et al. 



33| and Eggels 



3J] , it has only more recently found 



app 



28 



ications for flows 



35 



36 



37 



38|. 



in different conflgurations and physical conditions - see for e.g. Refs. 

The effects of particle collisions in the solution of LBE are generally represented by 
relaxation-type models. One of the most common among them is the single-relaxation-time 
(SRT) model, also termed as the Bhatnagar-Gross-Krook (BGK) model [3^. Owing to its 
simplicity, the use of the SRT model in the LBE has been popular for simulating a 

variety of problems, including the computation of turbulent flow problems mentioned above. 
It is well known that the SRT model is quite susceptible_to numerical instabilities when it is 
employed for simulating high Reynolds number flows [40i]. In particular, the lack of proper 
mechanisms to properly dissipate unphysical small-scale oscillations arising due to non- 
hydrodynamic or kinetic modes in the LBE can often cause numerical instabilities 4l|. In 
the case of turbulent flows, and more speciflcally in coarser grid eddy-capturing simulations, 
such spurious oscillations may interfere with turbulent fluctuations and can result in loss of 
accuracy and stability. An important approach to enhance numerical stability with using 
SRT models is through the entropic lattice Boltzmann methods (ELBM), which ensures 



positivity of the distribution functions 42 



43 



44 



45| . While being endowed with elegant 



and desirable physical features, it may be noted that they have certain computational and 



46, 



471], and are not the pursued in this current 



physical limitations, as pointed in Refs. 
work. 

A more general form of the LBE, sometimes also called the moment method or the gener- 



alized lattice Boltzmann equation (G 
(MRT) to represent collision effects 



jBE), is based on the use of multiple relaxation times 
481 ]. It is actually a reflned form of the quasi- linear 



cuaiiy a r( 

BQQ, 



relaxation version of LBE with a collision matrix [6|, |7|, |49|, where collision is carried out 
in the moment space. In contrast to the SRT-LBE, the MRT-LBE or GLBE deals with 
moments of the distribution functions, such as momentum and viscous stress directly. This 
moment representation provides a natural and convenient way to express various relaxation 
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processes due to collisions, which often occur at different time scales. Also, the collision 
matrix takes a much reduced form as a diagonal matrix in this moment space. By carefully 
choosing and separating different time scales to represent changes in various hydrodynamic 



and Mnetic „,odes through a von Neun,ann stabUity a„a,y.s of the Mnetic equation ^, 

the numerical stability of the LBE can be significantly improved [40|. The general forms 
of the MRT models in two-dimensions (2D) and three-dimensions (3D) are presented by 



5l[ |. respectively. Simplified forms of MRT 



Lallemand and Luo [40| and d'Humieres et al. 

models 52,[53,[5^ and with a different weighted representation of moments [55|| have also 
been introduced to improve boundary conditions and to improve ability to represent hydro- 
dynamics with thermal fluctuations. The MRT-LBE has been further extended with the use 
of additional forcing terms to simulate complex fluid flows, such as multiphase flows in 2D 
and 3D by McCracken and Abraham 56| and Premnath and Abraham 57|] , respectively, and 



applied to simulate complex multiphase flow problems with significantly enhanced numerical 



stabihty [58|, l59|. More recently, the GLBE approach [571] has also been used to simulate 



complex magnetohydrodynamic problems with much success 601 ]. 

In recent years, Yu et al. [61] developed a MRT-LBE for LES of certain classes of turbulent 
flows. In particular, they employed the Smagorinsky SGS model with a constant coefficient, 
where the local strain rate tensor is given in terms of non-equilibrium moments. As such, 
their approach is applicable for problems without boundary effects, for e.g., free-shear flows 
and they have indeed validated it for a turbulent free-jet flow problem. However, the presence 
of either stationary or moving boundaries, such as walls or free surfaces, respectively, are 
known to strongly affect the turbulence structures, and suitable modiflcations are needed to 
the standard Smagorinsky SGS model for use with the GLBE. Moreover, in many situations, 
external forces, such as constant body forces mimicking pressure gradient in a periodic 
domain or non-uniform forces such as Lorentz or Coriolis forces, can drive and/or strongly 
influence the character of turbulent flow physics. The effects of these forces can be introduced 
as forcing terms in the GLBE. Also, the use of forcing terms representing non-uniform forces 
provide a framework to introduce more general forms of SGS Reynolds stress models that 
are not based on eddy-viscosity assumption. Moreover, as the scales of turbulent flow vary 
locally in general situations, it is important to employ local grid reflnement approaches in 
conjunction with the MRT-LBE. 

Thus, a primary objective of this paper is to develop a framework for LES using the MRT- 
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LBE with forcing term for wall-bounded flows, in which near-wall turbulence is generally 
known to be anisotropic and inhomogeneous in nature. We propose to carry out the forcing 
term in the natural moment space of the GLBE so that it is readily amenable for simulating 
general forms of non-uniform forces. The computations of the moment-projections of the 
forcing term are provided for the three-dimensional, nineteen velocity (D3Q19) model ^. 
To account for the reduction in the turbulent length scale near walls, we employ the van 
Driest wall damping function 6^ in the Smagorinsky SGS model. We derive expressions 
for the strain rate tensor used in the SGS model in terms of the non-equilibrium moments 
of the GLBE in the presence of forcing terms representing general non-uniform forces by 
means of Chapman-Enskog analysis 57|, |63|], which is a generalization of those presented 
by Yu et al. [61]. We also briefly discuss an optimized computational procedure for such an 
extended GLBE formulation. Moreover, we incorporate variable resolutions in the GLBE by 



introducing a conservative local grid refinement approach 6J, |65||. While the use of a con- 
stant Smagorinsky SGS model is known to have certain limitations (see, for e.g., Ref. [sil), as 
a first step to model bounded flows as well as for reasons of computational efficiency, we have 
employed it in conjunction with the damping function, which are known to be reasonably 



accurate for certain wall-bounded flows 



66|. It may be noted that more sophisticated SGS 



models involving the use of dynamic procedures to determine the values of the parameters 
in the SGS mo dels H t.- ccu^veuts some of the limitations of the constant Smagorin- 
sky SGS model have also been successfully used in the LBM context recently by Premnath 



et al. 



68(1 . Another important recent development is an inertial-consistent Smagorinsky SGS 



model proposed for use with the LBM by Dong et al. 

Another objective of this paper is to perform systematic studies for assessment of accu- 
racy and gains in numerical stability using the LES framework described above for a set 
of canonical wall-generated flow turbulence problems. In particular, we evaluate the GLBE 
in detail for two problems viz., fully-developed turbulent channel flow at a shear or friction 
Reynolds number Re of 183.6 based on channel half-height and 3D driven cavity flows at 
a Re of 12, 000 based on cavity side width. The benchmark problems involve complex fea- 
tures of wall-bounded turbulent flows, and extensive prior data, including those from DNS 
and experiments are available to compare and assess the results of the detailed structure of 
turbulence statistics obtained using the GLBE computations. We also study the gains in nu- 
merical stability when the GLBE is used in lieu of the SRT-LBE for such complex anisotropic 
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and inhomogeneous turbulent flows as well as the parallel scalability of its implementation 
on a massively parallel cluster. 

This paper is organized as follows. In Section [ITl we discuss the development of the 
generalized lattice-Boltzmann equation (GLBE) with forcing term. Section UTTl presents the 
subgrid scale model for wall-bounded turbulent flows used in this work. Details of the com- 
putational procedure of GLBE and its optimization are provided in Sec. IIVI The simulation 
results, accuracy and stability of two canonical problems, viz., fully-developed turbulent 
channel flow and 3D cubical cavity flow are discussed in Sees. |V] and IVIt respectively. Fi- 
nally, the summary and conclusions are presented in Sec. IVIIi More elaborate details of the 
approach used in this work are presented in various appendices. 



II. GENERALIZED LATTICE BOLTZMANN EQUATION WITH FORCING 
TERM 



The computational approach for turbulent flows based on the solution of the GLBE is a 
recent version of the LBM. The GLBE consists of the evolution equation of the distribution 



function of particle populations as they move and collide on a lattice and is given by 



51 



57| 



faC^ + e^St,t + 6t) - fa{'x,t) = - ^A„^ (^f^- f^l"^ +^ Ma/3 - ^^a/?) ^ff^t- (1) 

13 /3 ^ ^ 

Here, the left hand side of Eq. ( [T]) corresponds to the change in the distribution function 
during a time interval St, as particle populations stream from location "af to their adjacent 
location ~x + e^St, with a velocity along the characteristic direction a. We consider a 
three-dimensional, nineteen velocity (D3Q19) particle velocities set, shown in Fig. [H given 

by 

(0,0,0) a = 

(±1, 0, 0), (0, ±1, 0), (0, 0, ±1) a = 1, ■ ■ ■ , 6 (2) 

(±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1) a = 7, ■ ■ ■ , 18. 

The magnitude of the Cartesian component c of the particle velocity is given by c = bxjSt-, 
where bt is the lattice time step. 

The first term on the right hand side (RHS) of Eq. ([1]) represents the cumulative effect 
of particle collisions on the evolution of distribution function /„. Collision is a relaxation 




FIG. 1: Schematic illustration of the three-dimensional, nineteen velocity (D3Q19) model. 

process in which relaxes to its local equilibrium value at a rate determined by the 
relaxation time matrix The GLBE has a generalized collision matrix with multiple 

relaxation times corresponding to the underlying physics: the macroscopic fields such as 
density, momentum and stress tensors are given as various kinetic moments of the distribu- 
tion function. For example, collision does not alter the densities p and momentum j = p'u , 
while the stress tensors relax during coUisins at rates determined by fiuid properties such 
as viscosities. The components of the collision matrix K^p in the GLBE are developed to 
refiect the underlying physics of collision as a relaxation process. 

The second term on the RHS of Eq.([l]) introduces changes in the evolution of distribution 
function due to external force fields F, such as driving body forces mimicking a pressure 
gradient in a periodic domain or gravity, or Lorentz or Coriolis forces, through a source term 
Sa- In this term is the component of the identity matrix /. The source term may be 
written as 0, Q 

S.= ^^^^^^^fT'\p,l!), (3) 



where f^^'^ip, 'u) is the local Maxwellian 



rj^^\p^it) = u^p{i + j^ + 



\ a = 



18 



a = l,---,6 (4) 



^ a = 7, ■ ■ ■ , 18. 



and Cg = cj \/3 is the speed of sound of the model. By neglecting terms of the order of 



O(Ma^) or higher Eq.([3]) may be simphfied as 



3 9 

Sa = Wa {Cai - Ui) + — {c^ ■ ll) tai Fi (5) 

where = {F;j., Fy, F^}, with F^, Fy and F^ are the Cartesian components of the external 
force field, which can, in general, vary in space and/or time. 

It may be noted that, Eq. ([T]) is obtained from the second-order trapezoidal dis- 



571], viz., + e^5t,t + 5t) - fa{'x',t) 



cretization of the source term in GLBE 
-Z]/3^a/3 [fi3{'x',t) - + l/2[Sa{'x',t) + Sai'x' + e^6t,t + 6t^]6t, which is made ef- 



fectively time-exphcit through a transformation fa = fa — ^f^SaSt [TOj, and then dropping 
the 'overbar' in subsequent representations for convenience. The second-order discretiza- 
tion provides a more accurate treatment of source terms, particulary in correctly recovering 
general forms of external non-uniform forces in the continuum limit without spurious terms 
due to discrete lattice effects [7l|, and its time-exphcit representation facilitates numerical 



solution in a manner analogous to the standard LBE. The local macroscopic density and 
velocity fields are given by 

a 

j ^ pit = ^fo,K + ^F6t, 

a 

and the pressure field p may be written as 



(7) 



P = CsP- 



(8) 



The physics behind the kinetic equation Eq. ([T]), and in particular, the collision matrix 
Aq,/3 will become more transparent when it is specified directly in terms of a set of linearly 
independent moments f instead of the distribution functions f = [/o, /i, • • • , /l8]^ i-e. 



through f 



/o, /i, /2, • • • , /: 



t 



18 . Here, the superscript 'f is the transpose operator and the 
'hat' represents quantities in moment space. The moments have direct physical import to 
the macroscopic quantities such as momentum and viscous stress tensor. The components 
of f are provided in Appendix [XI This is achieved through a transformation matrix T: 
f = Tf . The elements of T are given in d'Humieres et al. 
orthogonal to every other row. The essential principle for its construction is based on the 
observation that the collision matrix A becomes a diagonal matrix A through A = TAT^^ 
in a suitable orthogonal basis which can be obtained as combinations of monomials of the 



51| . Each row of this matrix is 



Cartesian components of the particle velocity directions through the standard Gram- 
Schmidt procedure. 

The collision matrix in moment space A may thus be written as 

A = diag {so, si, S2, . . . , Sis) , (9) 

where Sq, Si, S2, . . . , Sig relaxation time rates for the respective moments. The corre- 
sponding components of the local equilibrium distributions in moment space f'^'^ = 
/q'', fi'^, /l*^, . . . , are functions of the conserved moments, viz., local density and mo- 
mentum fields, and are given in Appendix |X] . 

When there is an external force field, as in Eq. represented in particle velocity space 
S, where S = [Sq, Si, 5*2, ... , •S'lg]^, appropriate source terms in moment space S need to be 
introduced. In this regard, we obtain the projections of source terms onto moments S by 
a direct application of the transformation matrix to Eq. ([5]) for the D3Q19 model, as done 

5*0, Si, S2, . . . , 5'i8 . 



for the D3Q15 model earlier in Ref. [57|], i.e. S = TS, where S 
They are explicit functions of the external force field F and the velocity Tt, which are 
summarized in Appendix El In effect, due to collisions and the presence of external forces, 
the distribution functions in moment space or simply, the moments are modified by the 
quantity -A (^f - f^«) + ^ ^/^^) 

A multiscale analysis based on the Chapman-Enskog expansion 63| of the GLBE shows 
that in the continuum limit, it corresponds to the weakly compressible Navier-Stokes equa- 
tions with external forces, where the density, velocity and pressure, given by Eqs. (El), (j7l) 
and (IH]), respectively, as was done for the D3Q15 model by Premnath and Abraham [57|. 
The macrodynamical equations can also be derived through an asymptotic analysis under 
a diffusive scaling 14, 0, 3- The transport properties of the fluid flow, such as bulk ( 
and shear u kinematic viscosities can be related to the appropriate relaxation times through 
either Chapman-Enskog analysis of the GLBE or the von Neumann stability analysis of its 
linearized version [4^: 

C ^ I - (10) 



9 Vsi 2 

^^ = ^(^-05i, /3 = 9, 11, 13, 14, 15. (11) 

Notice that from Eq. (fTTl) . sg = su = S13 = S14 = S15 to maintain isotropy of the stress tensor 
and S2 determines the magnitude of bulk viscosity. The rest of the relaxation parameters do 

10 



not affect hydrodynamics but can be chosen in such a way to enhance numerical stabihty as 
to simulate higher Reynolds number problems for a given grid resolution, in particular for 



ity analysis [40|], the 
5l|: si = 1.19, S2 = 



wall-bounded turbulent flows considered here. Based on linear stabi 
following values for the other relaxation parameters are determined 
Sio = = 1.4, S4 = Sg = Sg = 1.2 and Sig = S17 = Sig = 1.98. For the conserved moments, 
the values of the relaxation parameters are immaterial as their corresponding equilibrium 
distribution is set to the value of the respective moments itself. However, with forcing terms 



it is important that they be non-zero 



56 



57|. For simphcity, we set Sq = S3 = S5 = S7 = 1.0. 



It may be noted that all relaxation parameters have the following bound < Sq, < 2. In this 
paper, we employ the above values for the relaxation parameters. Since the GLBE employs 
a set of relaxation times, it is also referred to as the multiple-relaxation time (MRT)-LBE. 

It may be noted that when all the relaxation times are equal, i.e., si = S2 = . . . = 
Sis = 1/r, where r is the relaxation time, the GLBE reduces to the single- relaxation time 
(SRT)-LBE ^, based on the Bhatnagar, Gross and Krook model 39|. Its popularity 
and appeal lies in its apparent simplicity. However, the GLBE has marked advantages 
when compared with the SRT-LBE: for a given resolution, the GLBE is significantly more 
stable numerically and more accurate for problems with anisotropy, with an insignificant 
additional computational overhead, thereby allowing access to a greater range of problems, 
particularly at higher Reynolds numbers, to be reached than possible with the SRT-LBE. 
This is demonstrated later for two problems involving wall-generated turbulent flows. 



III. 



SUBGRID SCALE TURBULENCE MODEL 



In this paper, we have incorporated the subgrid scale (SGS) effects in the GLBE through 
the standard Smagorinsky model [32]. It assumes that the SGS Reynolds stress term depends 
on the local strain rate tensor and leads to the eddy- viscosity assumption. The eddy viscosity 

arising from this model can be written as 



(12) 



where Cg is a constant (taken equal to 0.12 in this paper). Here, A is the cut-off length 
scale set equal to the lattice-grid spacing, i.e. A = 5x, and Sij is the strain rate tensor given 
by Sij = 1/2 [djUi + diUj). In LBM, the strain rate tensor can be computed directly from 
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the non-equilibrium part of the moments, without the need to apply finite differencing to 



the velocity field. Recently, Yu et al. 



6l[ | derived such expressions for the strain-rate tensor 



for the D3Q19 model of the GLBE without forcing term. In this paper, we extend the 



results for GLBE with forcing term by means of a Chapman-Enskog analysis [57|, |63|, which 
is presented in Appendix [Bl The use of forcing terms allows for incorporation of not only 
general forms of non-uniform external forces, but also more general forms of SGS Reynolds 



stress models 



681]. This procedure for calculation of strain rates in GLBE is fully local in 



space and is computationally efficient, particularly for complex geometries. 

The eddy viscosity ut is added to the molecular viscosity z/q, obtained from the statement 
of the problem, through the characteristic dimensionless number, such as shear Reynolds 
number for turbulent channel fiow problem, to yield the total viscosity v (i.e., z/ = z/q + z/^). 
The relaxation times may then be obtained from Eq. ffTTl) . When such eddy- viscosity 
type SGS models are used to provide additional contributions to the relaxation times, the 
GLBE can be considered to be "coarse-grained" and it can be readily shown that the macro- 
scopic dynamical equations of fiuid fiow corresponds to the filtered equations with the SGS 
Reynolds stress represented through the eddy viscosity. As a result, the GLBE would repre- 
sent the dynamics of larger eddies in turbulent fiows. The distribution functions (or equiv- 
alently, the moments) and the hydrodynamic fields, can be considered to be grid-filtered 
quantities. An alternative approach is to directly apply filters to the moment representa- 
tion of the GLBE and rigorously derive SGS models essentially from kinetic theory under 
appropriate scaling 13]. ^ ^ 

To account for the damping of scales near the walls, following an earlier work ||66i], we 
have implemented the van Driest damping function 62 ] 



A = 5, 



1 — exp ( — 



v4^ 



(13) 



where is the distance from the wall and is a constant, taken to be 25 in this work [66]. 
The superscript + signifies normalization with respect to wall units, i.e. z'^ = z/5y, where 
5u = u^/vq is the characteristic viscous length scale. Here, u^, is the shear or friction velocity. 



which is related to the wall shear stress through = a/ t^^/ p. While this approach has 
some empiricism built-in, for a class of wall-bounded turbulent fiows, such as turbulent 



channel fiows considered here, it has been shown to be reasonab 
based on the solution of grid-filtered Navier-Stokes equations 



y accurate in prior work 



661 ]. Also, as will be shown 



12 



later in this paper that the GLBE is able to reproduce turbulence statistics in the near- 
wall region reasonable well using this damping supplemented to the SGS model. For more 
general situations, for improved accuracy it may be necessary to introduce dynamic SGS 



models (e.g., |67|, [T^, l75|) for LES using the GLBE |68|. 



IV. COMPUTATIONAL PROCEDURE: OPTIMIZATION WITH FORCING 
TERMS 



In practice, implementation of the GLBE with forcing term, i.e. Eq. ([T]), together with as- 
sociated turbulence models and procedure for strain rate computations, initial and boundary 
conditions, requires careful consideration for the details for efficient performance. In particu- 
lar, the "effective" collision step including the forcing terms should be performed in moment 
space, while the streaming step should be executed in particle velocity space and the special 



properties of the trans: 
fully exploited 



51 



57 



'ormation matrix that transform between the two spaces should be 



6l| . Such properties of the transformation matrix T include its or- 



thogonality, entries with many zero elements, and entries with many common elements that 
are integers, which are used to form the most compact common sub-expressions for trans- 
formations between spaces. We will now briefly discuss the details of the computational 
procedure for the GLBE with forcing term used in this paper. 

The GLBE with forcing term can be re-written in terms of the following "effective" 
collision and streaming steps, respectively: 



f (af, t) = f ( a?, t) + -cu( af , t). 



and 



where fa is the post-collision distribution function and 



t37( = 



-A f - f 



eq 



X--A1S 



(14) 



(15) 



(16) 



is the effective change due to collision including the effect of external forces. Here, f = 
f(af,t), P'^ = f^^(~?,t) and S = S(af,t), and X is the identity matrix and A = TAT^^ = 
diag{so, si, . . . , sis) is the diagonal collision matrix in moment space. 
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A note regarding the actual implementation details is in order. First, the transformation 
matrix T is row-wise orthogonal and satisfies TT''" = £, where T"^ is the transpose of T and 
£ is a diagonal normalization matrix. Thus it follows that the matrix inverse is obtained 
simply using T^^ = T'^ C As a result, we may write Eq. flTBl) as -a7("af,t) = T'''q where q 
is given by q = —V — f^^^ + {^^^ ~ S and F = £^^A. Thus, for computational 
efficiency, we actually implement the "effective" collision step that also including forcing 
terms in moment space. Now, the relaxation times in A used to compute in Eq. f|T6|) 
can be related to the transport coefficients and modulated by eddy viscosity, in the case 
of hydrodynamic time scales, as follows: s^^^ = s^^ = 1^ + 1 from Eq. ( |T0|) . and sg = 
s\\ = si3 = si4 = Si5 = s^, where = Su + ^ = 3(z/o + Ut) + ^, from Eq. ([TT]). The 
eddy viscosity ut is obtained from Eq. ( |T2l) . The rest of the relaxation parameters can 
be chosen to enhance numerical stability, as discussed in Section [Tll The forcing term 
used in the computation of strain rate tensor (Appendix [B]) and in the "effective" collision 
step (Eq. f fT6l) ) can be obtained from Appendix [Al This optimized procedure dramatically 
improves the computational speed of the GLBE as compared to a naive implementation. 
Indeed, the additional computational overhead of using GLBE in lieu of the SRT-LBE is 
small, between 15%— 30%, but, as will be shown later, with a significantly improved accuracy 
and numerical stability. 

No slip wall boundary conditions, involving stationary walls as well as moving walls, in 
the case of turbulent channel flow and driven cavity flow, respectively, are implemented by 



means of the link or half-way bounce back [52] • To initiate turbulence, a three-dimensional 
perturbation velocity field satisfying divergence free condition 
of the GLBE through the consistent initialization procedure 



^ is employed in the solution 

3- 



V. FULLY-DEVELOPED TURBULENT CHANNEL FLOW 

First, we simulated a canonical problem, viz., fully-developed turbulent channel fiow 
using the GLBE with the SGS model mentioned above. Prior efforts have validated LBM 
as a DNS tool for this problem by comparing a set of turbulent statistics with available 



data 
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25 



30| . The focus of this study is to evaluate MRT-LBE that incorporate subgrid 
scale effects for this problem on a relatively coarse grid, while maintaining the necessary 
near-wall resolution. 
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Slip surface 




FIG. 2: Schematic of computational domain for LES of fully-developed turbulent channel flow. 

We considered turbulent flow with a shear Reynolds number Re,, = u^^H/vq = 183.6, 
where z/q is the molecular kinematic viscosity and H is the channel half height. A schematic 
of the problem set up is shown in Fig. [2], in which a no-slip boundary is imposed at the 
bottom, free-slip at the top and periodic boundary conditions were applied in the streamwise 



X and spanwise y directions 76|, 1781]. The computational domain is chosen with appropriate 



aspect ratios, viz., 6H and 3H in the streamwise and spanwise directions, resp ectively. 
With this domain, a sufficient number of wall-layer streaks are accommodated 76| and end 
effects of two-point correlations are excluded, i.e. the two-point velocity correlations in 
solutions are required to decay nearly to zero within half the domain [tqI. For this initial 
case, we considered a uniform grid with a grid spacing in wall units (referred to with a "+" 
superscript) as A+ = A/5j, = 4.08, where 6i, = z/q/u* is the viscous length scale as defined 
in Sec. IIIII The computational domain thus consists of 270 x 135 x 47 grid nodes. Due 
to the use of link-bounce back method for implementation of wall boundary condition, the 
first lattice node is located at a distance of = A+/2, which in our case is 2.04. For 
wall-bounded turbulent flows, it is important to adequately resolve the near-wall, small-scale 
turbulent structures, which is satisfied when the computations resolve the local dissipative or 
Kolmogorov length scale r] = (z/q/c)^/^, i.e. A+^ < 0(?7+) [3]. In particular, it is generally 
recognized that 1.5?7+ — 2.0?7+ represents the upper limit of grid-spacing, above which the 
small scale turbulent motions in bounded flows are not well resolved. It can be shown by 
simple arguments that ri~^ ^ 1.5 — 2.0 at the wall and that ■i]'^ increases with increasing 
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distance from the wall 



17| . Thus, our computational set up is expected to fairly resolve the 



small-scale turbulent structures. 

The initial mean velocity is specified to satisfy the 1/7*^ power law 17|], while initial 
perturbations satisfying divergence free velocity field The density field is taken to be 
P = Po = 1-0 for the entire domain. The precise form of the initial fields may not affect 
affect the turbulence statistics, but can have significant influence on the number of time 
steps needed for convergence of the solution to statistically steady state. In particular, the 
above choice of initial fields would enable a rapid convergence to the statistically steady 
state solution of the fluctuating fields obtained by the GLBE with forcing term. With 
these initial conditions on the macroscopic fields, we employed the consistent initialization 
procedure for the distribution functions or moments 1?^]. Using F = —^x = = ^x 
as the driving force, the GLBE computations are carried out until stationary turbulence 
statistics are obtained, as measured by the invariant Reynolds stresses profiles. This initial 
run was carried out for a duration of 50T*, where T* = H/u^ is the characteristic time scale. 
The averaging of various flow quantities was carried out in time as well as in space in the 
homogeneous directions, i.e. over the horizontal planes, by an additional run for a period of 
SOT*. 

Figure [3] shows the computed mean velocity profile, normalized by the shear velocity u^,, 
as a function of the distance from the wall given in wall units, i.e. = z/S^,, where 6^, is the 
viscous length scale defined earlier. Also plotted are the DNS data by Kim, Moin and Moser 
(1987) jsO] based on a spectral method and the von Karman log-law of the wall, which 
is valid for the so-called log-region. The computed velocity profile follows the DNS data 
fairly closely, with about 5% difference. Such differences are characteristic of LES, which 
employ relatively coarser grids than DNS, and they also generally depend on the numerical 



dissipation of the computational approach for LES (see e.g., Ref. [8ll. |82[|). 

Figure H] shows comparison of the computed mean velocity profile as a function of the 
distance from the wall with wall-layer scaling laws, i.e. viscous sublayer and log-law of 
the wall. Generally, Reynolds stress effects are negligible in the viscous sublayer region 
{z~^ < 5), and =< u > /u^ = holds for the mean velocity. For z'^ > 30, the mean 
velocity satisfies the log-law, i.e. = A\nz~^ + B, where the coefficients depend on the 
flow parameters and nature of the wall. Values of A = 2.5 and B = 5.5 are known to be 



reasonably accurate for flow over smooth walls at Re* ^180 IT], 



76 



801]. It can be seen that 
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FIG. 3: Comparison of the computed mean velocity profile with wall-layer scaling laws in outer 
wall coordinates for fully-developed turbulent channel flow at Re^, = 183.6. 

computations agree well with these scaling laws. 

The Reynolds stress, normalized by the wall-shear stress, is presented in Fig.[5]in semi-log 



scale and compared with the DNS data of Kim, Moin and Moser [80|] obtained from the direct 
solution of incompressible Navier-Stokes equations (NSE), which the GLBE computations 
reproduce quite well with good accuracy. 

Let us now consider the statistics of turbulent fluctuations of important quantities in 
the near-wall region. Figures El [7] and M show comparisons of the components of the root- 
mean-square (rms) streamwise, spanwise and wall-normal velocity fluctuations, respectively, 
computed using the GLBE with data from DNS based on the solution of NSE by Kim, Moin 
and Moser 



801 1 and experimental measurements of Kreplin and Eckelmann j83|. It may be 
seen that the computed results agree reasonably well with prior data. 

Another important quantity representing turbulent activity near the wall is the pressure 
fluctuations. Figure [9] shows the computed rms pressure fluctuations. The proflle shown 
here is qualitatively consistent with the NSE-DNS results. It is found that the pressure 
fluctuations normalized by the wall shear stress is about 1.66 at the wall, which is within 
the range in prior data - NSE based DNS results in Ref. js^ and [3] provide values of about 
1.5 and 2.15 respectively. These values depend on the Reynolds number employed. In the 



measurements reported by Willmarth 8J] , the values of maximum rms pressure fluctuations 



17 




FIG. 4: Comparison of the computed mean velocity profile with wall-layer scaling laws in inner 
wall coordinates for fully-developed turbulent channel flow at Re,, = 183.6. 

were found to be between 2 and 3, but these were for much higher Reynolds numbers than 
considered here. Moreover, the computed maximum pressure fluctuations occurs at ~ 26, 



which is close to the range in DNS data 80[, i.e. z'^ ^ 30. It may be seen that there the 



computed profile rms pressure fluctuations using GLBE is systematically somewhat larger 
than the DNS results based on NSE. This observation is also consistent with those found 
in Ref. [3^, where DNS using SRT-LBE revealed similar values for the peak rms pressure 
fluctuations and its location. Such difference could plausibly be due to compressibility effects 



inherent in LBM, while the DNS carried out in Ref. [8OI] considered incompressible NSE. 

A particularly stringent test is the comparison of computed components of near wall rms 
vorticity fluctuations with DNS, which is shown in Fig. (TUl in which lines represent the 



GLBE solution and symbols the DNS data 



801 ]. The components of vorticity fluctuations 



normalized by the mean wall shear (uI/uq). There is a strong variation among the com- 
ponents of vorticity due to inhomogeneity and anisotropy of turbulence closer to the wall. 
Also, as expected, for distances further from the wall, all the components of vorticity are es- 
sentially the same. While there are some deviations from the DNS results, the GLBE is able 
to reproduce the qualitative trends and, more importantly, the ratio between components 
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FIG. 5: Reynolds stress normalized by the wall shear stress for fully-developed turbulent channel 
flow at Re* = 183.6. 



of vorticity. Such deviations have been observed in LES using filtered NSE, e.g. [85|, who 
point out the underestimation of the resolved components to be due to mere consequence of 
filtering process inherent to LES. 

Another important measure is the pressure-strain (PS) correlations. Their components 
are: PS^ = {pd^u^), PSy = {pdyu'y), and PSz = (pdzU^), where the prime denotes 
fluctuations and the brackets refer to averaging (along homogeneous spatial directions and 
time). They provide indications of energy transfers among the components. Figure [TT] shows 
the components of the computed PS correlations. They exhibit the expected behavior close 
to the wall, including the transfer of energy from the wall-normal component to the other 
two components near the wall - a phenomenon termed as splatting or impingement 66|. 
Thus, it appears that the GLBE with forcing term is a reliable approach for computation 
of fully-developed turbulent channel flows. 
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FIG. 6: Root-mean-square (rms) streamwise velocity fluctuations normalized by the wall shear 
velocity for fully-developed turbulent channel flow at Re* = 183.6. 
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FIG. 7: Root-mean-squarc (rms) spanwise velocity fluctuations normalized by the wall shear ve- 
locity for fully-developed turbulent channel flow at Re* = 183.6. 
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FIG. 8: Root-mean-square (rms) wall-normal velocity fluctuations normalized by the wall shear 
velocity for fully-developed turbulent channel flow at Re* = 183.6. 
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FIG. 9: Root-mean-squarc (rms) pressure fluctuations normalized by the wall shear stress for 
fully-developed turbulent channel flow at Re* = 183.6. 
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FIG. 10: Components of root-mean-square (rms) vorticity fluctuations normalized by a factor given 
in terms of wall shear stress for fully-developed turbulent channel flow at Re* = 183.6. Lines are 



GLBE solution and symbols are DNS data 
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FIG. 11: Components of computed pressure-strain correlations normalized by a factor given in 
terms of wall shear velocity and molecular viscosity for fully-developed turbulent channel flow at 
Re* = 183.6. 
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A. Numerical Stability 



To put things in perspective, let us now discuss the stabihty characteristics of the GLBE 
in relation to the SRT-LBE for turbulent channel flow on coarser grids. The test case used a 
shear Reynolds number Re,, of 180, with a uniform spacing A"*" of 6 in wall units, which at the 
near-wall node becomes A^^ = 3 due to the link-bounce back scheme employed. The number 
of grid nodes used in each case is 180 x 90 x 32. This is a somewhat coarser resolution than 
used in the previous simulation and it is expected that small-scale near-wall dynamics may 
not be properly resolved. Nevertheless, subgrid scale motions are quite energetic for such 
coarse resolutions and it is important to determine if the grid scale numerical instabilities 
developed by the computational approaches interact with them. The numerical stability of 
the LBM depends on various factors including the grid resolution A, maximum velocity or 
Mach number Ma considered and the relaxation times or the molecular viscosity of the fluid 
uq. For a given resolution and maximum flow velocity, the numerical stability of the LBE 
depends mainly on the molecular viscosity of the fluid z/q. 

As is natural for LBE, unless otherwise specified, all the results are reported in lattice 
units. That is, the velocities are scaled by the particle velocity c and the distance by the 
streaming distance of the populations, 6x. Here, we considered a maximum velocity, i.e. 
velocity at the top surface to be about 0.18, and varied the viscosity z/q. In the case of the 
SRT-LBE, the only parameter that can be used to specify is the single-relaxation time 
r and its value is chosen from uq = l/3(r — 1/2). On the other hand, for the GLBE, the 
relaxation parameters that determine moments involving fluid stresses are determined from 
Eq. [m while the rest of the parameters are tuned to improve numerical stability as specified 
earlier. 

Figure [12] shows the components of rms turbulent fluctuations obtained by using both the 
GLBE and the SRT-LBE at z/q = 0.0012. The rms turbulent fluctuations results from the 
SRT-LBE simulation show some physically unrealistic behavior, with a large spike in the wall 
normal component near the no-slip wall. Farther out, ripples which grow as the slip-surface 
is approached can be seen in both the wall normal and the streamwise component. That is, 
spurious oscillations due to non-hydro dynamic or kinetic modes seem to strongly interact 
with fluctuating turbulent motions generated by the wall, particularly in the wall-normal 
component, in the case of the SRT-LBE. In contrast, due to scale separation of relaxation 
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FIG. 12: Comparison of the components of root-mean-square (rms) velocity fluctuations normalized 
by the wall shear velocity for fully-developed turbulent channel flow with a free-slip surface at the 
top for Re* = 180 obtained by GLBE or MRT-LBE (dashed lines with open symbols) and BGK- 
LBE or SRT-LBE (solid line with filled symbols) on a coarse grid. 

parameters in the GLBE, the kinetic modes are quickly damped and do not exhibit such 
unphysical behavior. 

It may be noted that such spurious effects do not seem to manifest with the SRT-LBE, 
when fine enough resolution is employed, as was also noticed, for e.g. in Ref. 30|. On 
the other hand, for the same resolution, if the viscosity is lowered further, the SRT-LBE 
becomes unstable. Stable and physically realistic solutions can be obtained only for viscosity 
greater than 0.0018 in this particular case. On the other hand, the GLBE seems to predict 
correct physical and smoother behavior for all the components of velocity fluctuations for 
viscosity of 0.0012 shown in Fig. [12] and up to 0.0006 in our work. For this specific problem 



we thus obtain enhancement in stability by a 
the observations made for other problems 51 



actor o: 
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about 3, which is consistent with 



61 



86|. Thus, it appears that 



the GLBE is superior in terms of both physical fidelity and stability on coarser grid LES 
simulations of anisotropic and inhomogeneous turbulent flows, when it is used in lieu of the 
SRT-LBE. We will also discuss more on the stability aspects when we discuss about the 
other canonical problem considered in this paper. 
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FIG. 13: Schematic of conservative local refinement using multiblock grids with GLBE. 
B. Conservative Multiblock Approach for Local Grid Refinement 

Close to a wall, length scales are very small, requiring a fine grid to adequately resolve 
turbulent structures. Use of a grid fine enough to resolve the wall regions throughout the 
domain can entail significant computational cost, and this can be mitigated by introducing 
coarser grids farther from the wall, where turbulent length scales are larger. One approach 
is to consider using continuously varying grid resolutions, using a interpolated-supplemented 



LBM 



87| that effectively decouples particle velocity space represented by the lattice and 

1 known that interpolation could introduce sig- 
401 ]. which could severely affect the accuracy of 



the computational grid. However, it is we 
nificant numerical dissipation, see for e.g. 



solutions involving turbulent fluctuations, as was confirmed in numerical experiments during 



the course of this work. Thus, we consider loca 



ly embedded grid refinement approaches, and 



in particular their conservative versions 64, l65|] that enforce mass and momentum conserva- 
tion. Similar zonal embedded approaches have been successfully employed in computational 



approaches based on the solution of filtered NSE for LES of turbulent flows 

Figure [13] shows a schematic of such a multiblock approach in which a fine cubic lattice 
grid is used close to the wall and a coarser one, again cubic in shape, farther out. In order 
to facilitate the exchange of information at the interface between the grids, the spacing of 
the nodes changes by an integer factor, in this case two. As well as using different grid 
sizes, the two regions use different time steps (time step being proportional to grid size), 
and the computational cost required per unit volume is thus reduced by a factor of 16 in the 
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coarse grid. Figure [T3] shows a staggered grid arrangement, in which nodes on the fine and 
coarse sides of the interface are arranged in a manner that facihtates the imposition of mass 
and momentum conservation. Different blocks communicate with each other through the 
Coalesce and the Explode steps, in addition to the standard stream- and- collide procedure. 
The details are provided in Chen et al. 6J] and Rohde et al. 65j, and here we very briefly 
present the essential elements in what follows. The Coalesce procedure involves summing 
the particle populations on the flne nodes to provide new incoming particle populations 
for the corresponding coarse nodes. Similarly, the Explode step involves redistributing the 
populations on the coarse node to the surrounding flne nodes. These grid-communicating 
steps used in the multiblock approach presented in Chen et al. 6^ were incorporated in the 
GLBE framework in this work. 

We performed fully-developed turbulent channel flow at the same shear Reynolds number 
as before, i.e. 183.6, with different blocks, viz., flne block near the wall and coarse block in 
the bulk bounded by top free-slip surface. For the flne grid, we used a resolution A^^^^ = 4 
in wall-units (with A,;[^ = 2 due to link-bounce back) and a resolution of A+^^j,g = 8 in wall 
units for the coarse grid. We used 256 x 128 x 17 grids for the flne block and 128 x 64 x 17 
for the coarse block, which corresponds to similar aspect ratios to that used in the earlier 
simulations. The initial run and averaging times used were similar to that for the uniform 
grid case, viz., SOT* and 30T*, respectively. 

The mean velocity and Reynolds stress proflles computed using the GLBE with locally 
reflned multiblock grids are compared with uniform grid solution, again computed using 
GLBE, along with the DNS data [8^ , in Figs. [TH and [15], respectively. Generally, good 
agreement between various simulations can be seen. Some differences between the DNS data 
and the LES results based on the GLBE noticed in these flgures are similar to those found in 
LES based on flltered NSE. The components of the rms velocity fluctuations in streamwise, 
spanwise and wall-normal directions are presented in Figs. [161 [13 and [TBI respectively. 
Again, the multiblock GLBE based LES results are fairly in good agreement with the uniform 
grid GLBE as well as the DNS data for various components of velocity fluctuations. It is 
found that the velocity fluctuations and Reynolds stress are somewhat sensitive to numerical 
artifacts arising near grid-transition regions, i.e. at the interface between flne and coarse grid 



blocks, where they are slightly damped. Similar features have been noted in Ref. 



65] when 



the multiblock approach is employed for computation of certain classes of flows, having 
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FIG. 14: Mean velocity normalized by the wall shear velocity for fully-developed turbulent channel 
flow at Re=K = 183.6. Lines are GLBE results on locally refined multiblock (broken) and uniform 



(solid) grids, and symbols are Kim, Moin and Moser's DNS data (1987) 



8C]. 




FIG. 15: Reynolds stress normalized by the wall shear stress for fully-developed turbulent channel 



flow at Re=K = 183.6. Lines are GLBE results on locally refined multib. 
(solid) grids, and symbols Kim, Moin and Moser's DNS data (1987) 



ock (broken) and uniform 
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FIG. 16: Root-mean-square (rms) streamwise velocity fluctuations normalized by the wall shear 
velocity for fully-developed turbulent channel flow at Re* = 183.6. Lines are GLBE results on 
locally refined multiblock (broken) and uniform (solid) grids, and symbols are Kim, Moin and 
Moser's DNS data (1987) [80;]. 

flow components normal to grid interfaces. On the other hand, when the rms pressure 
fluctuations computed using the multiblock GLBE are compared with uniform grid results, 
there is a slight overprediction by the former, plausibly due to added compressibility effects 
with the use of multiblock grids (see Fig. [T9l) . 



C. Parallel Scalability 

One of the main advantages of the LBM is its natural amenability for implementation on 
parallel computers. The code implementation of GLBE with forcing term was parallelized 
using the Message Passing Interface library through a domain decomposition strategy that 
exploits the local and explicit nature of the approach. It was tested for parallel scalability 
on a large parallel cluster known as Seaborg located at U.S. Department of Energy's NERSC 
center. In these tests, the size of each subdomain was held constant at 20 x 256 x 256 (1.3 
million grid nodes), per processor for wall-bounded turbulence simulations. The speed-up 
factors obtained for up to 1024 processors are shown in Fig. [201 It is evident that near-linear 
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FIG. 17: Root-mean-square (rms) spanwise velocity fluctuations normalized by the wall shear 
velocity for fully-developed turbulent channel flow at Re* = 183.6. Lines are GLBE results on 
locally refined multiblock (broken) and uniform (solid) grids, and symbols are Kim, Moin and 
Moser 's DNS data (1987) [80,]. 

scaling can be obtained on massively parallel clusters and thus it appears that the GLBE 
approach is well-suited for large-scale turbulent flow simulations. 



VI. THREE-DIMENSIONAL FLOW IN A CUBICAL CAVITY 



Let us now consider another wall-bounded flow problem, viz., 3D flow in a cubical cavity 
driven by its top lid, and its computation using the GLBE. Although the geometry is simple, 
it is characterized by richness in fluid flow physics as there are no homogeneous directions 
and the presence of walls on all sides profoundly modifies the flow behavior. Features such as 
multiple counter-rotating recirculating regions at the corners, Taylor-Gortler-type vortices, 
bifurcations in flows and transition to turbulence may manifest themselves depending on the 



Reynolds number 



89|. Generally, when the Reynolds number Re based on the cavity side 



length is less than 2000, the flow field is laminar, and flow instabilities manifest themselves 
near the downstream corner eddy when Re is between 2000 and 3000. As Re increases, 
turbulence is generated near the cavity walls, with the flow near the downstream corner eddy 
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FIG. 18: Root-mean-square (rms) wall normal velocity fluctuations normalized by the wall shear 
velocity for fully-developed turbulent channel flow at Re* = 183.6. Lines are GLBE results on 
locally refined multiblock (broken) and uniform (solid) grids, and symbols are Kim, Moin and 
Moser's DNS data (1987) 



becoming fully turbulent when Re>10,000. Due to various states exhibited by the flow at 
higher Re it is a very challenging problem to study, in particular in obtaining computational 
results as it requires accurate methods with long averaging times. Measurements from 
experiments in cubical cavity are available for Re = 10, 000 in Prasad and Koseff 9^, while 
pseudo-spectral DNS and spectral element LES were performed more recently at Re 
12, 000 by Leriche and Gavril 



of LBM, d'Humieres et al. 



akis 9l| and Bouffanais et al. 93] , respectively. In the context 
5l[ | performed simulations of 3D flow in a diagonally driven cavity, 
in the laminar and transition regime, i.e. Re < 4000. The focus here is to perform GLBE 
simulations at a higher Re of 12, 000 and compare results with available prior computational 
results 



91 



921 ] and experimental data [90] . In addition, we also compare numerical stability 



of the GLBE and the SRT-LBE for this problem at higher Re range. 

The computational conditions that we considered are as follows. The schematic of the 
3D flow in a cubic cavity of side length 2W is shown in Fig. [2T] with a coordinate system, 
in which flow is driven by the top lid with velocity Uq. The Reynolds number used in 
our computations Re = Uq2W/ z/q = 12, 000 is achieved by setting the lid velocity to be 
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FIG. 19: Root-mean-square (rms) pressure fluctuations normalized by the wall shear stress for 
fully-developed turbulent channel flow at Re* = 183.6. Lines are GLBE results on locally refined 
multiblock (broken) and uniform (solid) grids, and symbols are Kim, Moin and Moser's DNS data 



(1987) 



Uq = 0.12, with a viscosity of uq = 0.00128 on a relatively fine uniform grid with 128'^ lattice 
nodes. A note regarding the choice of the lid velocity is in order. For a given Re, when the 
fluid viscosity is chosen based on relaxation parameter so as to maintain numerical stability, 
the choice of Uq influences the number of grid nodes needed to resolve each side of the 3D 
cubic cavity. That is, any reduction of Uq by a factor k will increase the total number of grid 
nodes by /c^. On the other hand, since the GLBE is a weakly compressible computational 
approach, the Mach number Ma (= Uq/cs where Cg = should be small. Thus, the 

value of Uq is chosen as a compromise between satisfying the weakly compressible condition 
and resolution requirements so as the obtain an acceptable level of solution accuracy. In this 
work, as found later, the choice of Uq = 0.12 and a resolution with 128^ grid nodes yields 
reasonably good accuracy. 

Now, imposing a constant lid velocity profile on the top lid leads to edge and corner singu- 
larities and can significantly affect the stability, convergence and accuracy of simulations at 



such high Re 



91 



92l | . In reality, there is a velocity distribution at the top lid, whose precise 



form is not known. Following Leriche and Gavrilakis QlJ as well as Bouffanais et at 92| . 
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FIG. 20: Parallel scalability of GLBE for turbulence simulation on a parallel cluster, Seaborg, 
at Department of Energy's NERSC. Symbols are GLBE performance and line is the ideal linear 
speed-up. 



we set the following velocity profile for the lid: 



uiid{x,y) = uo 
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It was found that the flow field in the cavity is not overly sensitive to the lid velocity profile, 



when such higher-order polynomial distributions are used [91|, |92| . The mean value of this 
velocity profile is Um ~ 0.85f/o, with over 75% area of the lid has a velocity above Um and 
the corresponding Reynolds number on the mean velocity is 10, 200. In the GLBE, the 
velocity boundary condition at the lid is provided by setting the distribution function of 
incoming populations corresponding to ~ea through an momentum-augmented bounce back 
as follows [93|: 



fa = fa + 2w„Po 



^Cf ' ^ lid 



(18) 



where ~ea = a- No-slip zero velocity boundary conditions based on bounce back ap- 
proach are set for all the other walls. A statistically stationary state of the flow field is 
obtained after running for SOOT* and then collecting statistics at each grid node averaged 
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FIG. 21: Schematic arrangement of computational domain for LES of flow in a three-dimensional 
cubical cavity of half- width W driven by its top lid with velocity Uq. 

over a period of 150T*, where the characteristic time is T* = 2W/uo- 

Figures [22] and [23] show the computed first-order statistics, viz., the mean velocity profile 
on two of the cavity centerlines along with the other available data for comparison. It is 



seen that GLBE solution is in reasonable agreement with the DNS 9l| and experimental 
data Q. 



In general, as discussed in Ref. [9l|, momentum transfer from the lid creates a region of 
high pressure in the upper corner of the downstream wall as the flow has to change direction, 
dissipating part of its energy. The flow then convects downwards along the downstream wall 
like an unsteady wall jet, which separates from the wall near the mid-section of the wall 
and leading to two elliptical jets. They subsequently impinge on the bottom cavity wall and 
generate turbulence, which is convected away by the central and main vortex. The second- 
order statistics of the fluctuating flow fleld provide an indication of the turbulent activity. 
Figures [21] and [25] provide the rms velocity fluctuations along the direction parallel to the 
lid motion, i.e. Urms on the centerlines y = W, z = W and x = W,y = W, respectively. 
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FIG. 22: Mean velocity < u > on the center line x = W,y = W obtained in LES using GLBE (solid 
line) compared with DNS of Leriche and Gavrilakis (2000) [91] (dashed line) and experimental 
data of Prasad and Koseff (1989) [90^ (circles). 
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FIG. 23: Mean velocity < w > on the center line y = W, z = W obtained in LES using GLBE (solid 
line) compared with DNS of Leriche and Gavrilakis (2000) 9lj] (dashed line) and experimental data 



of Prasad and Koseff (1989) 



90] (circles). 
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FIG. 24: Root mean square (rms) velocity fluctuations urms on the centerline y = W, z = W 
obtained in LES using GLBE (solid line) compared with DNS of Leriche and Gavrilakis (2000) 



(dashed hne), LES with NSE of Bouffanais et al. (2007) 



ml 



92l ] based on dynamic model (dotted 



line) and dynamic mixed model (dot-dashed line) and experimental data of Prasad and KosefF 



(1989) 



90(] (circles). 



while Figs. [21] and [27] provide the rms velocity fluctuations along direction normal to the lid 
motion i.e. Wrms on the centerlines y = W, z = W and x = W,y = W, respectively. In 
addition, the components of the Reynolds stress < u'w' > on the centerlines y = W, z = W 
and X = W,y = W are provided in Figs. [28] and [291, respectively. Firstly, these results 
indeed show that turbulence is generated along cavity walls. In particular, the turbulent 
fluctuations are about an order of magnitude larger near the downstream wall than near 
the upstream wall. Moreover, the fluctuations are the largest along the bottom wall. These 
seem to be consistent with the description of the features of the fluid motion in the cavity, as 



elucidated in the DNS 



91| . Although there is some deviation in the peaks of the fluctuations 



S and LES considered here are based 



91 



92l |. the GLBE computations, in 



when compared with other data, considering that DNS 
on approaches using higher-order spectral methods 
general, compare reasonably well with them, which are very encouraging. 

Some differences observed between computed solutions, including those from DNS and 



LES 



91 



92l |. and the experimental data could be attributed to differences in Reynolds 
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FIG. 25: Root mean square (rms) velocity fluctuations urms on the centerline x = W,y = W 
obtained in LES using GLBE (solid line) compared with DNS of Leriche and Gavrilakis (2000) 



ml 



92l ] based on dynamic model (dotted 



(dashed hne), LES with NSE of Bouffanais et al. (2007) 
line) and dynamic mixed model (dot-dashed line) and experimental data of Prasad and Koseff 



(1989) 



90(] (circles). 



number as well as the averaging times used. For example, the magnitude of the peak value 
of the near-wall Reynolds stress in Fig. [28] is influenced by the length of the time interval 
over which averaging is performed. In this work, we have chosen the time period of averaging 
(150T*) from the sampling period used in experiments 90|]. In effect, our results are closer 
to these data. However, it should be noted that prior computations [oil, [o^ found that the 
peak value of Reynolds stress is conditioned by rare events, which occur on time intervals of 
approximately SOT*. Hence, the averaging period in Refs. 9ll,l92] is chosen such that the rare 



events, which tend to suppress fluctuations, are sampled many times, which accounts 



or the 



92|. 



difference between experiments |90|] (including this work) and prior computations |9lL 

Moreover, a note regarding the influence of the choice of the SGS turbulence model on 
the turbulence statistics results is in order. Bouffanais et al. 92], who employed dynamic 



SGS models in their computations yielding high fidelity results in excellent agreement with 
DNS, also reported preliminary results without employing a SGS model, i.e. unresolved 
DNS, and with a constant Smagorinsky SGS model. They found that the unresolved DNS 
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FIG. 26: Root mean square (rms) velocity fluctuations wrms on the centerline y = W,z = 
W obtained in LES using GLBE (solid line) compared with DNS of Leriche and Gavrilakis 



(2000) 



based on dynamic model 



91[ (dashed line), LES with NSE of Bouffanais et al. (2007) [93] 
(dotted line) and dynamic mixed model (dot-dashed line) and experimental data of Prasad and 
Koseff (1989) [9Q] (circles). 

is not even qualitatively correct for this problem and the use of constant Smagorinsky SGS 
model resulted in improved predictions, but still not fully consistent with the resolved DNS 
and their results with using dynamic models. However, unlike Bouffanais et al. [92], in 



this work, we have used van Driest wall damping function [62| in conjunction with the 
constant Smagorinsky SGS model, which appears to make considerable difference in further 
improving the quality of the results. It appears that the use of wall damping function, which 
accounts for reduction of near- wall turbulent length scales, yields results that are significantly 



closer and more consistent than without t 
consistent with our recent observation 



le use of such a damping function. This is also 
68| that the use of a wall damping function with 
constant Smagorinsky SGS model results in better agreement with the dynamic SGS model 
results for LES of turbulent channel flow than without it. It should, however, also be noted 
that while the results obtained with the use of damping function are generally better, they 
are in some cases somewhat underpredicted near walls, e.g. Fig. [251 ^ind the peaks in some 
other cases are broader, e.g. Fig. [261 While the motivation for the use of damping function 
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FIG. 27: Root mean square (rms) velocity fluctuations wrms on the centerline x = W,y = 
W obtained in LES using GLBE (solid line) compared with DNS of Leriche and Gavrilakis 

based on dynamic model 



(2000) loi] (dashed line), LES with NSE of Bouffanais et al. (2007) c 
(dotted line) and dynamic mixed model (dot-dashed line) and experimental data of Prasad and 
Koseff (1989) [90] (circles). 



is simplicity and computational 



models in the LBM framework [6^] for further improvements 



efficiency, it is, of course, desirable to employ dynamic SGS 



A. Numerical Stability 

We will now make direct comparisons of stability characteristics of the GLBE with the 
SRT-LBE for 3D cavity flow simulations at higher Reynolds number ranges, complementing 
an earlier study ^] . In both approaches, for a given grid resolution, the shear viscosity was 
fixed and the lid velocity was increased gradually until the computation became unstable. 
Figure [30] shows the maximum Reynolds number that could be attained before the compu- 
tations became unstable. Results are provided for different grid resolutions and viscosities 
for both the approaches. The superior stability characteristics of the GLBE are evident for 
this wall-bounded turbulent flow problem. The GLBE computations can reach Reynolds 
numbers that are several times higher than that of the SRT-LBE, typically by a factor of 
3 or 4 and sometimes even about as high as an order magnitude. Indeed, the SRT-LBE 
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FIG. 28: Reynolds stress < u'w' > on the center line y = W,z = 
(solid line) compared with DNS of Leriche and Gavrilakis (2000) 



V obtained in LES using GLBE 



of Bouffanais et al. (2007) 



9l|] (dashed line), LES with NSE 



921 ] based on dynamic model (dotted line) and dynamic mixed model 



(dot-dashed line) and experimental data of Prasad and Koseff (1989) [90] (circles). 



became unstable and unable to simulate the Re = 12, 000 case considered above with a 
resolution of 128^ using the GLBE. 
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FIG. 29: Reynolds stress < u'w' > on the center line x = W,y = 
(solid line) compared with DNS of Leriche and Gavrilakis (2000) 



V obtained in LES using GLBE 



of Bouffanais et al. (2007) 



91[ (dashed line), LES with NSE 



921 ] based on dynamic model (dotted line) and dynamic mixed model 



(dot-dashed line) and experimental data of Prasad and Koseff (1989) 90(] (circles) 
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FIG. 30: (Color online) Comparison of numerical stability characteristics of GLBE with SRT- 
LBE for 3D cavity flows: Maximum attainable Reynolds number for a given resolution and shear 
viscosity. 
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VII. SUMMARY AND CONCLUSIONS 



A generalized lattice Boltzmann equation (GLBE) with forcing term, which uses multiple 
relaxation times, for eddy-capturing computations of wall-bounded turbulent flows that are 
characterized by statistical anisotropy and inhomogeneity is discussed. Standard Smagorin- 
sky eddy viscosity model is used to represent SGS turbulence effects, which is modified by 
the van Driest damping function to account for reduction of turbulent length scales near 
walls. Second-order and effectively time-explicit source terms, which represent general forms 
of non-uniform external forces that drive or modulate the character of turbulent flow, are 
projected onto the natural moment space of GLBE in this formulation. In this framework, 
the strain tensor used in the SGS model is related to the non-equilibrium moments and 
the forcing terms in moment space. Furthermore, local grid refinement using a conservative 
multiblock approach is used to coarsen grids in the bulk flow region, where turbulent dissipa- 
tion or Kolmogorov length scales become larger. Computational optimization, particularly 
in the presence of moment-projections of the forcing terms, is also discussed. 

Two canonical bounded flows, viz., fully-developed turbulent channel flow and 3D driven 
cavity flow have been simulated using this approach for shear Reynolds number of 183.6 and 
Reynolds number based on cavity side length of 12, 000 respectively. The structure of turbu- 
lent flow given in terms of turbulence statistics, including mean velocity and components of 
root-mean-squarc (rms) velocity and vorticity fluctuations and Reynolds stress are in good 
agreement with prior DNS and experimental data. The computed rms pressure fluctuations 
are found to be somewhat over-predicted in comparison with DNS data, which is based on 
the solution of incompressible Navier-Stokes equations. It is thought this may due to the 
kinetic nature of the GLBE approach, which is inherently weakly compressible. In the case 
of 3D cavity flow, the GLBE is able to capture turbulent velocity fluctuations and Reynolds 
stresses generated by cavity walls, which are in reasonably good agreement with prior data. 

With regard to numerical stability, it is found that by separating various relaxation 
times, the GLBE is able to maintain solution fidehty, while the SRT-LBE solution can 
exhibit spurious effects on velocity fluctuations in the near-wall region, particularly in the 
wall normal component, on relatively coarser grids in turbulent channel flow simulations. 
The GLBE is found to be superior in maintaining numerical stability at higher Reynolds 
number Re for 3D cavity flows, with the maximum attainable Re several times that for 
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the SRT-LBE, depending on the resolution and shear viscosity of the fluid. Moreover, 
parallel implementation of the GLBE approach is able to maintain near-linear scalability in 
performance for over a thousand processors on a large parallel cluster. 

The GLBE with forcing term appears to be a reliable approach for LES of wall-bounded 
turbulent flows. It is expected that further improvements can be achieved by introducing 



n 



more advanced SGS models based on dynamic procedures [67|, 0, Q in the LBM 68|. 
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APPENDIX A: COMPONENTS OF MOMENTS, EQUILIBRIUM MOMENTS 
AND MOMENT-PROJECTIONS OF FORCING TERMS FOR THE D3Q19 LAT- 
TICE 



The components of the various elements in the moments are as follows 51 1: /o = P) /i = 

e,f2 = e^,f3 = j^, U = Qx, % = jy, h = Qy, fr = jz, fs = Qzjg = ^Pxx, fio = ^t^xx, fn = 

Pww, fl2 = T^ww, /l3 = Pxy, /l4 = Pyz, fib = Pxzi /l6 = ^xi /l7 = "f^yi fl8 = ""^z- 

Here, p is the density, e and represent kinetic energy that is independent of density 
and square of energy, respectively; jx, jy and jz are the components of the momentum, i.e. 
jx = pUx, jy = pUy, jz = pUz, Qx, %, Qz are the components of the energy flux, and Pxx, 
Pxy, Pyz and Pxz are the components of the symmetric traceless viscous stress tensor. The 
other two normal components of the viscous stress tensor, pyy and Pzz, can be constructed 
from Pxx and Pww, where p^w = Pyy — Pzz- Other moments include Hxx, t^ww, ^x, ^y and 
rriz- The flrst two of these moments have the same symmetry as the diagonal part of the 
traceless viscous tensor pij, while the last three vectors are parts of a third rank tensor, with 
the symmetry of jkPmn- 
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The components of the equihbrium moments for the D3Q19 lattice are as follows: 51| : 



/o^^ = p,/f^ ^ e^^ = -Up + 19^, ^ e^'^" = 3p - 'i^Jf = J.JT_^ = q7 = 

_2 ■ feq _ ■ Teq ^ eq _ _2 ■ feq _ ■ 7eq ^ eg _ _2 ■ ^eg ^ o eg _ [^-^'^ ~ ' ] fe? ^ 

o„eq _ o ( _l^eq\ feq _ eg _ [-^l"-^^] feg _ eg _ _l„eg Jeg _ eg _ jxjy_ 7eq _ eg _ 
"-"'xx ~ I, 2t'xx) 'ill — /^ujw ~ p yJl2 — ''ww ~ 2^wwyJl3 — Pxy ~ p ) i 14 — Pyz ~ 

'^Jlt^pll = tJt! = o,/;i = 0, A1 = 0. 

The components of the source terms in moment space can be obtained by multiplying 
the transformation matrix with Eq. ([5]). The final expressions are as follows: 

5o = 0, Si = 38{FxUa; + FyUy + FzUz), S2 = -lliF^u^ + FyUy + F^u^), 

S3 = Fx, 5*4 = — f-Fc, S'5 = Fy, Sq = —^Fy, Sy = Fz, Ss = ~^Fz, 

Sg = 2{2FxUx - FyUy - F-zU-z), Siq = -(2F^m^. - FyUy - FzUz), 

Sll = 2{FyUy - FzUz), S12 = -{FyUy - FzUz), S 13 = {F ^U y ^ FyU x) , 

Si4, = {FyUz + FzUy), Si5 = {FxUz + FzUx), S^ie = 0, Sn = 0, Sis = 0. 

The self-consistency of the moment projections of source terms is evident. For e.g., S3, S^ 
and 5*7 provide Cartesian components of body forces on the moments corresponding to the 
components of momentum (mass flux), Si provides the work due to forces on the moment 
corresponding to kinetic energy, etc. 



APPENDIX B: STRAIN RATE TENSOR USING NON-EQUILIBRIUM MO- 
MENTS IN THE GLBE WITH FORCING TERM 

In this section, we will present a brief derivation of the strain rate tensor in terms of the 
non-equilibrium moments of the GLBE with forcing term by applying a Chapman-Enskog 



analysis. The results that follow are generalizations of those presented by Yu et al. [6l|| to 
include forcing terms representing non- uniform forces. First, the left hand side of the GLBE 
is simplified by applying Taylor series, which results in the following: 



-T-'A (f - P^) + (^l - S6t, (Bl) 



where = diag (^dt^ , dt^ + et ■ V , . . . , dt^ + eis ■ . Now applying the Chapman-Enskog 
expansion 

f^f^^ + fW^t (B2) 
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to Eq. (IBip . and truncating terms of order 0{S^) or higher, we get 



= -Af (1) + ( I _ Ia ) s, (B3) 



where = T~BtT~^ = diag (^dtg,dtg + Eiidi, . . . ,dtg +Eisidi^, in which and henceforth 
summation of repeated indices is assumed. Eq. ( 1B3I) can be rewritten in terms of non- 
equihbrium moments f-^^ as 



f ~ -A-' [dt, + E,a) +(a- S, 



(B4) 



where E^i = Te^jT"^. 

Now substituting the expressions for the equihbrium moments P'' and the source terms S 
in Eq. flB4p . we simphfy the expressions for the components of the non-equihbrium moments. 
Some such components of interest are as follows: 

fl'' - e(^) = -i Id,, (-Up + 19^-^1 + Idu^ +(--l]s,, (B5) 



i2 a2 



- pS - ~ {9.,. (^) \ (4.. + ^'...)} H- - i) S., (B8) 

" ^ (^) + 5 <^»^-' + ^'^4 + 4) (ss) 

/iS' ^ p£> = -i- |a„ j + 1 (aj, + aj,,) I + (^_L _ (bio) 

For further simplification, we invoke the following approximations: dt, 
2n,.F,, dt, (I) ~ 2uyFy, dt, (f ) ~ 2u,F,, dt, (^) - u^Fy + dt, 
Mj^F^ + M^Fj^, dtQ - ^^x-F^ + u^F^, and 9t(,p = -Sfcjfc, which result in 

= -f - ^^1, (Bll) 

^ - - iSg, (B12) 
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J 11 


2 1 

3 sii 


{dyjy - dj^) - 




J 13 ^ 


1 1 

3 Si3 


(dxjy + dyj^) - 


2^13 


f(l) ^ 
il4 ^ 


1 1 

3 Si4 


(dyjz + djy) - 


-6l4 


ftl) ^ 
J 15 ^ 


1 1 

3 Si5 


{dxjz + dzjx) - 


2^15 



(B13) 
(B14) 
(B15) 
(B16) 

It follows from the above that the components of the strain rate tensor can be written 
explicitly in terms of non-equilibrium moments as 

1 



St. 



yy 



S. 
5. 



38p 
1 

1 

'76p 



'■11 



19 i^Sgk 

19 (sg/Jr" + 3sii/in / 



xy 



2p 
3 



(neq) 
13 ' 



n4 5 



2p 
2p 



'15 



where 



f eg J c 

2 



a G {1,9,11,13,14,15} 



(B17) 
(B18) 
(B19) 
(B20) 
(B21) 
(B22) 

(B23) 



Here, the components of the source term Sa can be obtained from App endix [Al The form 



of Sij turns out to be very similar to that obtained by Yu et al. [6l|, except for the ex- 
pression h^a^'^\ which contains the additional contribution |S'q, that provides the effect of 
the forcing term. The procedure discussed here, however, is general, and can be read- 
ily employed for deriving the expressions for strain rate tensor for other lattice velocity 
models in the presence of forcing terms. The magnitude of the strain rate IS"! used in 



turbulence models can then be obtained from Eqs. flBT7l) - flB22|) as |5| = v^25~5^ 



2(5'|a; + S'^y + + 2(5*2^ + 5*^2 + Sl^)). To clarify the notations employed, we again note 
that 5*0, represents the source terms in moment space, Sa corresponds to the relaxation times 
in the collision term, and Sij is the strain rate tensor. 
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